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Abstract 

Nonlinear mixed effects models represent a powerful tool to simultaneously analyze 
data from several individuals. In this study a compartmental model of leucine ki- 
netics is examined and extended with a stochastic differential equation to model 
non-steady state concentrations of free leucine in the plasma. Data obtained from 
tracer /tracee experiments for a group of healthy control individuals and a group of 
individuals suffering from diabetes mellitus type 2 are analyzed. We find that the 
interindividual variation of the model parameters is much smaller for the nonlinear 
mixed effects models, compared to traditional estimates obtained from each indi- 
vidual separately. Using the mixed effects approach, the population parameters are 
estimated well also when only half of the data are used for each individual. For a 
typical individual the amount of free leucine is predicted to vary with a standard 
deviation of 8.9% around a mean value during the experiment. Moreover, leucine 
degradation and protein uptake of leucine is smaller, proteolysis larger, and the 
amount of free leucine in the body is much larger for the diabetic individuals than 
the control individuals. In conclusion nonlinear mixed effects models offers improved 
estimates for model parameters in complex models based on tracer /tracee data and 
may be a suitable tool to reduce data sampling in clinical studies. 
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1 Introduction 



The powerful combination of nonlinear mixed effects (NLME) mode ls and stochas t ic dif- 



ferential equations (SDEs) h as lately received increasing attention ( Donnet et al. . 2010l : 



modelling (lOvergaard et all . l2005l : iM^Uer et all . |2010| ). 



Picchini and Ditlevsenl. 120101) , with many applications in the area of pharmacokinetic/pharmacodynamic 



Nonlinear mixed effects models are used to study population characteristics when data 
have been gathered for multiple i ndividuals governed by the same intraindividual (within 
individuals) mechanisms (see, e.g..lDayidian and Giltinanl.ll995U2003t lOlofsen et al.l . 12004 ; 
Pillai et ali . 120051 : IPinheiro and Batesl . I2OOOI : ISheiner and Beat Il980h . NLME models are 
hierarchical in structure and take population behaviour and individual properties into 
account simultaneously. They provide a statistical framework to separate intra- and in- 
terindividual (between individuals) variability by using probabi lity distributions f or the 
individual-specific parameters. This stands in contrast to what ISheiner and Beall (jl980|) 
call the two stage approach to population modelling, where model parameters are esti- 
mated separately for each individual. 

Therefore NLME models give better estimates of the variabilities and it has been shown 
that the interindividual variabilities of the parameters are better estimated for NLME 



models than for two stage models (see, e.g., lOlofsen et al.l . 12004 ISheiner and Beall . Il980 



19831 : ISteimer et al.l . Il984l ). Moreover, as long as the number of individuals is sufficiently 



large, a smaller amount of data samples is needed for each indiv idual in order to obtain 
good estimates of population properties. This was, e.g., shown by lOvergaard et al.l (|2005[ ) 
in a simulation stu dy of a one-cornpartm ent pharmacokinetics model based on a single 
SDE. Furthermore, iJonsson et al.l ( I2OOOI ) show that NLME models perform better than 
two stage models to detect and characterize nonlinearities in the physiological system. 

NLME models are traditionally based on a function that is nonli near in the population 
(fixed effects) parameters and in the individual-specific parameters (jPavidian and Giltinan 
I995I ). A term is also added to accommodate for noise in the measurements. Ordinary dif- 
ferential equations (ODEs) are often used to represent state variables in dynamic models. 
The model output may be a function of the state, time, model parameters, and a measure- 
ment noise term. To accommodate for noise that is intrinsic to the system, and not only 
due to measurement errors, a better approach may be to use stochastic differential equa- 
tions. Such intrinsic noise may be due to true randomness in the model parameters and 
state variables over time. Also in systems without intrinsic noise, stochastic differential 
equations can be used to reduce model errors due to an oversimplified ODE model. 

In this study compartmental models and tracer/tracee experiments are analyzed using 
techniques from nonlinear mixed effects modelling. As an example we focus on a model 
of leucine kinetics based on tracer/tracee data. Figure [1] shows a graphical representation 
of the model. In several publications this model, or minor modifications thereof, has been 
used as a part of larger models, in which the kinetics of lipoproteins in human blood 
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plasma has been studied (see, e.g., iBarrett et al.l . l2006l . for a review). However, to our 
knowledge, it is the first time the model is used in an NLME setting. 

The evolution over time of the model state variables is described by a system of four 
coupled ODEs. In this paper we also extend the ODE model with a mean-reverting 
Ornstein-Uhlenbeck process to model a non-steady state concentration of free leucine in 
the plasma through the experiment, which turns the system of ODEs into a system of 
stochastic differential equations, SDEs. As far as we know, it is the first time this approach 
is used on tracer/tracee data. Moreover we believe that there is no published material, 
up to date, where NLME and SDEs are combined in such a complex model (w.r.t. the 
number of model state variables and unknown parameters to estimate). 

Experimental data from 19 control individuals and 15 individuals suffering from diabetes 
mellitus type 2 (DM2) are used to estimate the parameters. Furthermore we investigate 
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larger ODE based compartmental model. 

We are aware of only a few studies that compare the t wo stage and NLME approaches 
for analysis of experimental data. An example is given in lOlofsen et al.l (120041 ). Therefore 
one main purpose of this paper is to compare different modelling approaches, and to show 
that nonlinear mixed effects models and the two stage approach give different estimates 
of population parameters. We are also interested in investigating the parameter estimates 
for the SDE model, and how much that model predicts the leucine level to fluctuate during 
the experiment. 




Figur e 1: A graphical repre sentation of the leuc ine kinetics model. It was developed bv lDemant et a I 
( 19961 ), based on the work in ICobelli et al.l ( 199ll ). Compartment 1 corresponds to the amount of free 
leucine in plasma. Compartments 3 and 4 correspond to a body protein pool with leucine uptake and 
slow release back to the plasma. Compartment 2 is a hepatic intracellular compartment from which 
the liver is fed with leucine that is used when apolipoproteins are synthesized. The fractional transfer 
coefficient between compartments j and i is denoted fcy , and Ui is the continuous inflow of material 
into the first compartment. A more detailed description of the model is given in Section [3TTI 
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1.1 Results in summary 



1.1.1 Comparison between the two stage approach and the NLME approach 



In our study the NLME models predict a more homogeneous population. This re- 
flects the known fact that the interindivi dual variability is o f ten overestimated when 



the two stage approach is used (see, e.g., lOlofsen et al.l . 120041 ISheiner and Beall . Il980 



19831 : ISteimer et al.l . Il984l ). We have also observed that, for some parameters, the two 



modelling approaches produced di fferent estimates o f popu lation averages. In simulation 
studies with large residual errors, ISheiner and Beall (119831 ) have shown that the NLME 
approach produce less biased estimates of population averages. 



1.1.2 Using a smaller data set for each individual 

Often rich sampling of data for each individual is difficult, time consuming, or expensive. 
One main advantage with NLME modelling is that information from multiple individuals 
is used simultaneously. Therefore, population parameters can often be estimated when 
the number of data points is small for each individual, which may be of importance in 
experimental planning. 

When only half of the data were used for each individual, the estimated parameters of 
the ODE based NLME model, where no noise inherent to the system is modelled for, 
were similar to when all data were used. For the SDE based NLME model the variability 
intrinsic to the system was estimated to be zero, indicating that data need to be sufficiently 
densely sampled when SDEs are used. For the two stage approach we could not estimate 
the parameters for the smaller data set. Thus, for the model investigated in this study, 
the data may be sampled less dense if the proposed ODE based NLME model is used in 
the analysis. 



1.1.3 Comparison between NLME models based on ODEs and SDEs 

When all parameters were allowed to vary between individuals both the NLME models 
based on ODEs and SDEs produced excellent fits to data and small correlations between 
residuals were observed. The estimated parameters were similar with the two approaches, 
but the interindividual variations were in general smaller for the SDE model. 

The amount of free leucine in the body fluctuates over time in the SDE model. When 
all individual leucine levels are normalized, by dividing with the estimated average for 
that individual, the standard deviation of the leucine level is estimated to be 0.089 at 
each time point. This means that, for a typical individual, the amount of free leucine is 
predicted to vary with a standard deviation of 8.9% of the estimated average value, which 
indicate daily variations of leucine levels. 



4 



1.1.4 Differences between populations 

For the two groups of individuals different estimations of the population medians of the 
parameters have been observed. Especially the level of free leucine was significantly higher 
for the diabetic subjects than the control individuals. Also the rate parameters differed 
between the two groups and the fraction of free leucine that was degraded or used for 
protein synthesis per time unit is smaller for the diabetic individuals. Moreover proteolysis 
is predicted to be larger for the diabetic individuals. 



2 Theory and methods 

In this section we use a simple one-state compartmental model to describe the theory 
of compartmental models, tracer /tracee studies, and nonlinear mixed effects models. We 
also describe how stochastic differential equations can be used to account for a fluctuating 
tracee in the system. Moreover we discuss how to define an objective function to evaluate 
how well a model output resembles data. 



2.1 Compartmental models - Tracer studies 

A compartment is de fined as a well mixed and kinetically homogeneous amount of material 



(jCobelli et al.l . 120001 ). A compartmental model describes the relation between the amount 
of matter in a compartment and the fluxes of material in and out from the system, but 
also fluxes between compartments (in case of multi-compartment models). 

Consider the one-compartment model described in Figure [21 The influx of material into 
the compartment per time unit is denoted U, the fraction of material that leaves the 
compartment per time unit is denoted a, and the amount of material in the compartment 
is denoted Q. A mathematical representation of the model can be expressed with an 
ordinary differential equation, 

dQ/dt = -aQ + U, (1) 
where, in the general, representation a and U may be state and/or time dependent. 



When the system is in steady state, i.e., when Q is a constant and = —aQ + U, 
measurements of Q will give little information about the dynamical properties of the 
system (unless a or U are known, it is then enough with only one measurement point of 



U 




Figure 2: A graphical representation of a one-compartment model. 
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Q to describe the complete system). A common way to gain more information about the 
system is by tracer experiments, in which a known amount of tracer material is added 
into the system. Thus, a new input signal, u = u{t), \s introduced for the tracer, g, and 
measurements of the tracer-to-tracee ratio, q/Q, at successive time points may be used to 
get estimates of a, [/, and Q. The tracer should h ave the same kinetic properties as th e 



tracee but not affect the overall system behaviour (lAndersonl . Il983l : ICobelli et al.l . |2000[ ) 



Under the assumption that Q is in steady state, the model takes the form of an ODE 
dg/dt = —aq + u{t) for the tracer, and the output (response) function, Y{t) = q(t)/Q. 
However, measurements are usually error-prone and it is common to add a term v = v{t) 
to the output function, accounting for measurement errors but also for model uncertainty 
and noise effects intrinsic to the system. Measurements are performed at discrete time 
points tk, k = 1, . . . ,d, and the v(tk):s are often assumed to be Gaussian (with expectation 
zero) and independent, i.e., v{tk) and v(ti) are independent for k ^ I. When the output 
differs significantly between different time points an alternative is to assume a proportional 
error model, e.g., a lognormal distribution Y{tk) = exp(f (t^)). Observe that, by 
taking logarithms on both sides, this may be transformed to an additive noise model. 



2.2 Modelling uncertainty and intrinsic system noise 

Stochastic differential equations may be used to account for noise effects that are inherent 
to the system and not due to measurement errors. Such effects may be caused by true 
random variations in the biological process. Moreover, in many cases the systems to be 
modelled are not fully understood, or too complex to be modelled exactly by systems 
of ODEs. The models may then be improved by adding a system noise term (turning 
the ODEs into SDEs) representing the lumped effect of all non-explicitly mechanistically 
modelled effects in the system. 

In a tracer/tracee model the tracee is often assumed to be constant. However, this is not 
exactly true in real biological systems through the time course of an experiment, but it 
may vary around a mean value Qq. This may, e.g., be modelled by the mean reverting 
Ornstein-Uhlenbeck process 

dQt = a{Qo - Qt)dt + aQodWt, (2) 

where Wt is the standard Wiener process, aQo determines the magnitude of the system 
noise (the reason to include Qo here is explained in the population model below), and a is 
the rate at which Qt reverts towards the mean Qo (given in fraction of material per time 
unit). In the start of the experiment we assume that Q is Gaussian with mean Qq and 
variance ((jQo)^/(2a), which is the asymptotic mean and variance of Qt regardless of the 
initial value. Then the process Qt is stationary and Gaussian with mean E{Qt) = Qq and 



covariance Cov{Qs,Qt) = (^gLe-"(*-«) for s < t f lKaratzas and Shrevel . [l99lh . Thus, at 
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each time point t, 



jcrQo) 
2a 



(3) 



If there is reason to suspect that the error terms, f (tfc), are not Gaussian and independent 
these terms may be spht into two parts, v{tk) = fi(tfc) + V2{tk), where V2 is assumed to 
be Gaussian and independent, representing the measurement errors, but vi accounts for 
noise intrinsic to the system and may be mod elled by a stochasti c differential equation. 



e.g., a zero-mean Ornstein-Uhlenbeck process (IM^Uer et al.l . |2010( ). 



Other approaches to account for system noise and uncertainty are to allow other parame- 
ters to vary randomly over time or to lump the noise and uncertainty into additive terms 
in the right hand side of the ODEs. Assuming these additive terms to be stochastic 
processes turns the ODEs into stochastic differential equations. 



2.3 Nonlinear mixed effects models with stochastic differential 
equations applied to tracer/tracee experiments 

When time series data come from a number of individuals in a population, the model pa- 
rameters are traditionally estimated separately for each individual in the studied popula- 
tion (independently of population knowledge). One can then compute population averages 
and other statistical properties from the parameters obtained for eac h individual. This 



appro ach to population modelling is called the "two stage approach" (ISheiner arid Beal 



m 

1980l ) or the "standard two stage approach" (jSheiner and Beall . ll983l : ISteimer et al.l . ll984| ) 



Another approach to model population behaviour is the statistical framework of nonlinear 
mixed effects modelling. NLME models are hierarchical in structure with a clear separa- 
tion of the intraindividual and interindividual variations. The interindividual differences 
are characterized by individual-specific values of the parameters, for which the distribu- 
tion of the population is assumed to be known. The intraindividual variation is defined 
as the individual measurement errors and system noise. 

In Equations (jl])-® we use the simple one-compartment model to illustrate how NLME 
models with SDEs may be constructed. It is a model where the tracee fluctuates around 
a mean steady value. The interindividual variations are seen in ((Tj) and ([H]), whereas the 
intraindividual variations are accounted for in Equations (jlj and ([6]). The population 
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model takes the form 



dQj = a{Ql - Qi)dt + aQlAWt, Q*(0) ^ N (Q^ ) , (4) 



2a 

dql = {-a^ql + u\t))dt, 5^(0) = 0, (5) 

Yl = \og{^)+vl (6) 

= aexp(r/^), (7) 
Q^ = Qoexp(ryy, (8) 

where z = 1, . . . , denotes the individual, subscript k = 1, . . . , denotes the variable 
value at sample time point for the z:th individual, and the measurement errors, vl, 
are assumed to be independent and Gaussian with expectation zero and variance Sk- 
The measurements have been log-transformed which results in an additive noise structure 
in (ED. 

The initial condition of the tracee state variable in Equation (jl]) is the asymptotic dis- 
tribution of Ql- If the tracer material is introduced into the system as a bolus input at 
time zero, an impulse function can be used for the tracer input function u\t). A mathe- 
matically equivalent model would be to use = and g*(0) = gg, where is the amount 
of introduced tracer material. The reason for using aQ^ as a diffusion constant in (jl]) is 
that the level of Ql differs between individuals. 

In ([7]) and ([8]) a* and Qq depend on the fixed effects parameters a and Qq, thought 
to represent a typical individual in the population, and the individual-specific random 
parameters r]^ = [vIjVqV- this paper the latter are assumed to be Gaussian over the 
population with expectation zero and covariance matrix Q = diag(a;^, cjq). Lognormal 
distributions are used for a* and Qq mainly since it prevents parameters from being 
negative, even if the interindividual variabilities are large. The parameters a and Qn do 



not r epresent means, but instead the medians of the lognormal distribution ([Blackwood 



19921). 



Note that th e model defined by (Ill)- (l8l) is similar to the pharmacokinetics NLME model 



with SDEs in lOvergaard et al.l (120051 ). There the plasma volume corresponds to the tracee 
variable, and is assumed to be constant. Noise is added to the differential equation that 
describes the kinetics of the drug, corresponding to the tracer Equation ([5]) in the above 
model. 



2.4 Approximation of the population likelihood function 

For parameter estimation, we apply the maximum likelihood framework where the goal is 
to find the parameter set that maximizes the likelihood function, defined as the probability 
of finding the measured output for a given set of parameters. In the presentation below 
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we allow for multiple measurement variables, in contrast to in where we only have one 
measured variable. 



To construct a likelihood function we assume that the densities in the output function 
are Gaussian at all time points, tk, and dependent only on the information available at 
time tk-i- For k = 2, . . . ,d\ let yi_i = [y{, ■ ■ ■ , vl^i]'^ be the data up to time t^^i for 
individual i. Let 

^EiY^\yi-i) (9) 



^ k\k~ 



be the predicted mean value at time tk, conditioned on the information available at time 
tk-i, and 



^i\k- 



Cov{Y^\yi 



k-lJy 



(10) 



yi 

^ k\k-l 



the corresponding conditional covariance matrix. Then the residual vector e 
is Gaussian with expectation zero and covariance i?^. The starting point, (^/|05 Ri), for the 
recursion algorithm needed to compute (pl)- (fTOl) is determined by the initial distributions 
of the state variables. 



Kalman filtering techniques may be used to compute and Rl 



The Kalman filter (KF) 
can be applied if the model is linear. It is optimal in the sense that it minimizes the vari- 
ance of the prediction error if it is assumed that the densities of the model outputs and 
the initial conditions are Gaussian, and that the covariance matrix for the Gaussian noise 



in th e state equations and the output noise covariance are state independent (iJazwinskil . 
1970l ). The basic idea is that at every sample the conditional distributions of the state 
variables are updated by taking the new data point into account. This typically results in 
a net movement from the predicted mean value of the state in the direction of the experi- 
mental observations. For non hnear models, such as the one in Equation (ED, the extended 



Kalm an filter may be used (IKristensen et al.l . l2004l : lOvergaard et al.l . l2005l : iJazwinski 



19701 ). It is based on a linearization (in each time step,tA;) of the nonlinear parts of the 
model, which results in approximate equations to which the KF can be applied. 

The population likelihood function is a marginal likelihood function, where all the indi- 
vidual parameters, ri\ are integrated out. The Gaussian densities of the residuals and 
the individual-specific parameters, 77*, gives a likelihood function of the form (see, e.g., 
Overgaard et al.ll2005l ) 



N 

m 3^) = n / ^MW, 



where 9 is the vector of population parameters (i.e., 9 = [a, Qo, ujI,Uq, a, S]'^ in model 

y = {y^, . . . , 3^^} is the complete set of measurement data. The individual log-likelihood 

function = /*(?]*) = l^{r]^; y\ 9) is given by 



1 d} 
k=l 



■ T 



e\ R- el + \og\Rl\+m\og{2ix) 



^r7' + log|l]| + rlog(27r) , (12) 
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where m is the number of outputs and r is the number of interindividual random effects, 
i.e., the dimension of rj^. 

In general it is not possible to find an exact expression of the population likelihood 
function for a model with SDEs. The integrand must therefore be approximated and we 
have applied the Laplace approximatio n method together with the first-or der conditional 
estimation (FOCE) method (see, e.g.. lOvergaard et al.ll2005t IWangl 120071 ). The Laplace 
approximation method is based on a second order truncation of the Taylor expansion of 
/* around a stationary point 17*. In the FOCE method 17* is estimated as 



?7* = argmin { -/*(■?]*; y\9)] 



(13) 



in co ntrast to the first order (FO) method that uses 17* = (see, e.g., lOvergaard et al. 



20051 . for the details). With the Laplace approximation method the population likelihood 
function in f lTT]) is approximated as 



N 



my) 



i=l 



f27r) 



-1/2 



(14) 



where, by disregarding second order derivatives of the residuals, e^, and covariance ma- 
trices, the Hessian Hr,{l\ff )) is approximated as 



driadr]b 



k=l 



de] 



dr], 



a orjb 




-Ki. t 



diib 



(15) 



k=l 



-Rl, 



I dRl 



-R 



-idRl 



dr]a drib 



in position {a,b), a,b = 1, . ..,r. Proving (1151) is eleme ntary by using differentiation 
rules for matrices found in [Petersen and PedersenI ( 120081 ). Differentiation with respect 
to a matrix (or vector) in (|T5|) is performed on all elements in the matrix. If there is no 
interaction between the interindividual effects, 77', and the intraindividual variation effects. 



the derivatives of Rl are zero and (if^(Z*(i7*)))^ ^ ^ 



drjb 



)a,b- 



Every separate evaluation of the approximate population likelihood function requires the 
estimation of the individual parameters r/* when the FOCE method is applied. This may 
result in computationally heavy calculations, especially when the number, r, of interindi- 
vidual effects is large. 
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2.5 Implementation 



All the implementations have been performed in MATLAB® (2008b, The MathWorks, 
Natick, Massachusetts, USA). 

2.5.1 Optimization methods 

For NLME models the population parameters are estimated as 

^ = argmin{-log(L(0;3^))}. (16) 
e 

To optimize ( TT^ we have used a quasi- Newton method ba sed on the so called BFGS (B roy- 



den, Fletcher, Goldfarb, and Shanno) updating formula (INocedal and Wrightl . Il999l ). To 



find the optimal population para meters in ([T6|) we sta r t with the derivative free Nelder- 



Mead simplex-reflection method (INocedal and Wrightl . Il999l ) for a finite number of iter- 



ations to get closer to the true optimum. Then the optimization has been fine-tuned by 
the BFGS-method. The gradients used in the BFGS method have been approximated by 
a finite central-difference approximation. 

The Nelder-Mead simplex-reflection method and the BFGS method are both local op- 
timization methods, i.e., only local minima of the objective functions are searched for. 
However, for the models investigated in this paper, different initial guesses of the param- 
eters always resulted in convergence to the same optimal parameter set, indicating that 
the objective functions do not have multiple local minima. 

The steps that we have used to estimate the parameters in the NLME settings are: 



• Choose an initial guess, 6*0, of the population parameters and rfQ = Q for the 
individual parameters (z = 1 . . . , A^). 

• Calculate the population likelihood ( IT^ by optimizing 77* in f ll3p for i = 1, . . . , A^. 

Let j7q denote the optimized rj^:s. 

• Set j = 1 and repeat the following until a stopping criterion is reached: 

- Find a new 6j according to the optimization algorithm by using as initial 
values in every new calculation of ( !T4|) . 

- Denote 17] the optimized r/* for 6 = 6j. 

- Set j=j + 1. 



11 



2.5.2 Inference on estimated parameters 



For inference on tlie estimated parameters of an NLME model an approximation of the 
covariance matrix may be obtained from 



Cov{9) ^ (Hei-logimy))))- 



(17) 



( iKristensen et al.l . 120041 ). where Hg denotes the Hessian with respect to 6. The maximum 
hkehhood estimate is Gaussian and the covariance matrix is the inverse of the Fisher in- 
formation matrix, which is the expected value of t he Hessian of th e negative log-likehhood 
function when data are seen as random variables f lPawitanl - lioOlh . For a given realization, 
3^, of the data, th e obse rved Fisher information matrix, Hg{— \og{L{9; y))), may be used 
instead (jPawitanl . 1200 ll ). However, in this study we only approximate the population 
likelihood function in (I14p and the Hessian is computed numerically by means of a central 
difference procedure around the estimated parameters. 

When different mathematical models are used it is interesting to test if they give sig- 
nificantly different estimated values of the same physical quantities. With the NLME 
approach we approximate the confidence intervals around each estimated parameter by 
using (fT7|) and assuming that 6 is Gaussian. Confidence intervals around the estimated 
population parameters can also be calculated by taking every individual into account, as 
discussed in Section 13.1.11 

If the confidence intervals around a parameter are disjoint when two different models 
are used, we say that the models predict the parameter differentially with statistical sig- 
nificance at a certain confidence level. In a study where sam ple means from rep eated 
measurements from normally distributed data were compared, iPayton et al.l (120001 ) sug- 
gest to use 85% confidence intervals to test the hypothesis that the means are equal at 
significance level 0.05. We follow that suggestion when we compare the parameters of 
our models, but it should be mentioned that the confidence intervals in our study are not 
calculated from repeated realizations from a normal distribution. 



2.5.3 The extended Kalman filter 



The extended Kalman filter has been implemented as described in lKristensen et al.l (j2004| ). 
The EKF needs estimates of the predicted initial state and the initial state covariance 
matrix. As discussed in Section fI7I\ we use Var{Q{0)) = {aQQy/{2a) as an estimate of 
the variance of the tracee state variable modelled as in Equation (jl]). When the underlying 
model is deterministic in the state variables (ODE based), the state covariance is zero. 

The state equations for the models examined in this study are linear differential equations 
and the only non-linear parts of the models are present in the model output equations. 
Therefore the linearization approximation in the EKF is only applied to the model output 
equation. 
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3 Population Models of Leucine Kinetics 



In this article NLME models based on ordinary and stochastic differential equations 
are constructed for a compartmental model that describes the kinetics of leucine in 
blood plasma. The ODE based model has been used in several publications as a sub- 



Adiels et al. 


( 


2005b 


); 



Barrett et al.l (120061 ) . As far as we know the model parameters are estimated on single 
individuals and population features have been inferred with the two stage approach in all 
studies up to date. 

When the leucine model has been applied to fit parameters to data from single individuals 
the estimated values of the parameters vary much over the population. An explanation 



for th is may be that measurement errors in data are not taken care of ( lOlofsen et al. 



20041 ). One benefit with the NLME approach is that it handles intraindividual variations 
explicitly. In our implementation large interindividual variations are penalized since the 
individual log-likelihood in Equation ( !T2|) becomes smaller for larger absolute values of 
the elements in r/* (through the term —r]'^^Vt~^r]^). 

A stochastic differential equation is used to account for tracee deviations from steady state 
as in Equation (jl]). It should be pointed out that a population approach such as NLME 
modelling is needed to estimate parameters in SDEs when the data are sparsely sampled 
for each individual, since rich sampling is needed t o separate systera noise from measure- 
ment noise in single subject estimation algorithms (lOvergaard et al.l . l2005l ). A population 
approach such as an NLME model allows for more data to be used simultaneously. 

The data used are taken from experiments where stable isotope-labelled leucine i^H^- 
leucine) was introduced in the system as a bolus input. The experimental setup is de- 
scribed in lAdiels et al.l (l2005al Jbl. l2006l ). The data set is divided into one group of control 
individuals and one group of individuals suffering from diab etes mellitus type 2. T he dat a 
from the 19 individuals in the control group are also used in lAdiels et al.l (2005a 0,|2pO6l). 
The data from the 15 diabetic individuals are also used in I Adiels et al.l (l2005al . 12006. ). 
Each individual was fasting 12 hours prior to the tracer injection and also during the 
sampling period. The samples were collected 16 times at 2, 4, 6, 8, 10, 12, 15, 20, 30, 
and 45 minutes and 1, 2, 3, 4, 6, and 8 hours after the tracer injection. However there are 
missing data for some individuals. The data for a single control individual is visualized 
in Figure 131 

We will now describe the ODE based individual model (additive Gaussian measurement 
noise) used in a two stage approach to indirectly compute population parameter properties 
(Section l3.ip . Then we present the ODE based NLME model (additive Gaussian mea- 
surement noise, lognormal distributed population parameters) used directly to estimate 
parameter values and their distributions taking a maximum likelihood approach (Sec- 
tion [3]2]). Thereafter the SDE based NLME model (additive system noise) is discussed 
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(Section 13.31) . The estimated population parameters are presented after each model has 
been presented. Finally we compare the results for the NLME models based on ODEs 
and SDKs (Section [321) • 



3.1 The ODE based individual model used in a two stage ap- 
proach 



Leucine is an essential amino acid, which means that it is not synthesized in the body 
and that humans must acquire it from the diet. It is an important element in apolipopro- 
tein B (apoB), and since low density lipoproteins contain exactly one apolipoprotein B-lOO 
molecule, stable isotope-labelled leucine can be used as tracer in kinetic studies of low den- 
sity lipoproteins. The model used in this article describes the kinetics of leucine in plasma 
before it enter s the a poB synthesis machinery in the liver. The model was developed by 
Demant et al.l ( 1l996l ). and is based on the work in lCobelli et al.l (1l99l[ ). The data consist 
only of measurements of free leucine in plasma and is expressed as a tracer-to-tracee mass 
ratio. Figure [1] shows a schematic view of the model. 

Compartment 1 corresponds to the amount of free leucine in the plasma. It has a con- 
tinuous inflow of leucine from other parts of the body and also an exit (corresponding to 
catabolism of leucine and uptake of leucine by body proteins with a very slow turnover 
rate), described by the fractional transfer coefficient parameter /cqi. The amount of inflow 
of unlabelled free leucine per time unit is given by the parameter Ui. Since the tracer was 
introduced as a bolus input all the labelled leucine is assumed to be in compartment 1 
at time zero. No additional labelled leucine is assumed to be introduced into the system. 
Compartments 3 and 4 correspond to a body protein pool that account for the uptake 
and subsequent slow release of leucine back to the plasma. Compartment 2 is a hepatic 



Data for individual 1 in tile controi group 




Figure 3: Logarithm of the data for an individual in the control group. The output predicted by the 
ODE based NLME model described in Section \32\ is also included in the figure. 
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intracellular compartment. From this compartment the liver is fed with leucine that is 
used when apolipoproteins are synthesized. The parameter that describes the fraction 
of material in compartment 2 that reach the apoB synthesis machinery in the liver per 
time unit is here denoted kL2- When time passes the system is washed out from labelled 
material and the amount of tracer in the system converges to zero. Throughout this paper 
the time unit used is hours, thus all fractional transfer coefficients are presented in the 
unit h~^. The amount of material in compartments is presented in mg. 



The o r dinary differential eq uations in the model of the tracee system are (lAdiels et al. 



2005bl : iDemant et al.l . Il996r ) 



dQi/dt = -{koi + hi + hi)Qi + h2Q2 + A^isQs + Ui (18) 

dQ2/dt = ~{ki2 + kL2)Q2 + k2lQl (19) 

dgs/dt = -{hz + ^43)Q3 + hiQi + A;34Q4 (20) 
dQJdt = -kuQi + k^aQs, (21) 

where Qj = Qj{t) is the amount of tracee material in compartment j at time t. We 
assume here that the fractional transfer coefficients are not state dependent. Hence the 
system is linear, and the systems for the tracer and the tracee are identical except for Ui 
which is controlled by the experimentalist for the tracer system (it is assumed that no 
labelled material is produced in the body). 

Let Q = [Qi, Q2, Qs, Q^]'^ and q = [qi, q2, qs, q^]'^ be state vectors for the tracee- and tracer 
system respectively, and U = [Ui, 0,0,0]^ be a vector of input variables in the tracee 
system. Then the complete model can be written compactly in matrix notation as 

dQ/dt = KQ + U (22) 
dq/dt = Kq (23) 
n = gi(tfc)/Qi(4), (24) 

where f l2^ corresponds to the tracee system, fl2^ to the tracer system, and ([21]) corre- 
sponds to model output taken at discrete time points t^, k = 1,. . . ,d, where d is the 
number of samples. Observe that the experiments are based on a bolus injection in com- 
partment 1 at time zero. It is assumed that the tracer material is instantly well stirred 
throughout the plasma and that there is no tracer material in the system before the ex- 
periments start, therefore ^(O) = [qio, 0, 0, 0]"^, where gio is the known amount of labelled 
leucine injected (7 mg/[kg body weight] in our case). Another assumption is that the 
exchange of hydrogen atoms between leucine molecules and other molecules in the body 
is negligible, i.e., it is assumed that labelled leucine molecules remain labelled through- 
out the experiment. The matrix K is the compartmental matrix of the system, i.e., the 
coefficient matrix of the linear system f lT8|) -f l2T|) . 



The tracee system is approximated by a steady state (lAdiels et al.l . l2005bl ). i.e., dQ/dt = 
and the fractional transfer coefficients, kij, are constant. In matrix notation we get 
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= KQ + U ^ U = —KQ or Q = —K~^U for the tracee equations. Compartmental ma- 
trices are always invertible if there are at least one exit compartment (a compartment with 
a flow out of the system) and the system contains no traps (subsystems with no flows to 
the outside of the subsystem) (ICobelli et al.l . |2000| ). 



When the leucine system is considered separately (without the apoB part) the fractional 
rate, (per hour), of material leaving compartment 2 for lipoprotein synthesis is uniden- 
tifiable and in this article it is set to O.Ol/i"^. This value has been chosen because 
it is the value that we have found for most subjects in an implementation of the full 
apoB-model. Further parameter constraints have been used to reduce the number of 
unknowns. This is actually necessary to get an a priori identifiable model since, with- 
out the const r aints, three different para meter sets des c ribe the same model structure 
(ICobelli et al.l . Il980l ). We have followed [Packard et al.l ( l2000l ) and used k2i = fci2 (the 



rate of the flow of leucine through the cell walls of liver cells are the same in both direc- 
tions) and = 0.1/^43. It is not necessary to estimate both Ui and Qi separately since 
U = -KQ =^ Ui = Qi{kQi{ki2 + kL2) + k2ikL2)/{ki2 + kL2) mg/hour. Thus, the parame- 
ters that are estimated for a single individual are /cqi, ki2, ki^, k^i, k^s, and Qi. 



3.1.1 Estimated parameters 

Assuming a steady state for the tracee, the parameters were estimated for each individual 
separately. We used the weighted least squares method, taking the measured values as 
weights at each time point. 

The geometric mean and the coefficients of variation (CoV) of the full population are 
presented in Tabled] for the six parameters. The reason for reporting the geometric mean 
instead of the arithmetic mean is that the parameters are approximately lognormally dis- 
tributed in the population. Therefore, the variations of the parameters are multiplicative. 
The geometr ic mean is a max imum likelihood estimator of the median of a lognormal 
distribution (jBlackwoodl . Il992l ) . The coefficient of variation of the multiplicative model is 
a normalized measure of the spread of the parameters. It relates the standard deviation to 
the median of the parameters and is calculated by CoV = 100 x A/exp(s2) — 1 %, where 
is the unbiased estimate of the variance of the logarithm of the estimated parameters 



f lBlackwoodl . ll992f ). 



The parameters ki2 and kis are the ones that display the largest variation over the popula- 





koi 


ki2 


ki3 


kai 


k43 


Qi 


All individuals 


2.15 (41%) 


1.77 (116%) 


1.78 (115%) 


3.48 (70%) 


1.00 (63%) 


373 (65%) 


Controls 


2.32 (32%) 


1.70 (137%) 


1.61 (132%) 


4.18 (37%) 


0.94 (53%) 


289 (43%) 


Diabetics 


1.96 (51%) 


1.86 (98%) 


2.04 (95%) 


2.76 (98%) 


1.09 (80%) 


516 (73%) 



Table 1: Estimated geometric mean (and coefFicient of variation) of the parameters when data from 
the 34 individuals are used. The population parameters are trimmed as explained in the text. 
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tion. This is, at least partially, because the objective function is less sensitive to variations 
in these parameters. Therefore they are more influenced by noise in the observations than 
the other parameters. The estimated population parameters are trimmed in the sense that 
a few estimates are manually removed from the calculations (e.g., fcoi is estimated to be 
0.001 - the predefined lower bound - for two individuals). All parameters except Qi had 
a few individual estimates that were manually removed. 

To compare the results from the two stage approach with the results from the NLME ap- 
proach in the forthcoming sections we calculated 85% confidence intervals for the medians 
of the parameters for respective group (see Figure HI where the corresponding intervals 
for the NLME models are also included). We especially observe that, for the two stage 
estimates, Qi differs between control group 1 and the diabetic group. The confidence 
intervals have been calculated as explained below. 

If a parameter p is normally distributed over the population with unknown mean and 
variance, a 100(1 — a)% confidence interval of the population mean of p can be estimated 
asp ± ta/2,N-i^^i where p is the sample mean of the N observations of p, s the (unbiased) 
sample standard deviation of the observations, and t^ /?. ;v_i is the 100(1 — a/2) percentile 



of the t-distribution with N — 1 degrees of freedom ( Larsen and Marx . 200ll ). However 



since the parameters are assumed to be lognormally distributed over the population we 
calculated confidence intervals (a, h) for the means of the logarithms of the estimated pa- 
rameters and u sed (e", e^) to obta in confidence intervals of the medians of the parameters 



as suggested in iBlackwoodI ( 1l992l ) . 



To compare the interindividual variations for the two stage approach and the NLME mod- 
els in Sections 13.21 and 13.31 we have also computed confidence intervals of the population 
variance of log(p*), where for an individual i is its actual parameter estimate divided 
by the group geometric mean of p for that individual. The confidence interval around 
o"^, can be calculated since is y^-di stributed with N — 1 d egrees of 



freedom if the sample comes from a normal distribution (ILarsen and Marxl . l200ll ). We as- 
sume that log(p*), i = 1, . . . , A^, is normally distributed and the estimated 85% confidence 
interval for the interindividual variations are presented in Figure |H 



3.2 The ODE based NLME model 

Equations ( I22|) -(l2^ and the above discussion lead to the following population model, 

dgi/dt = ir^gl, g^ = feo,0,0,Or, (25) 

r,^ = log(^)+.;^, vl^N{Q,S), (26) 

where superscript i = 1, . . . ,N denotes the r.th individual. When considering one individ- 
ual, this is the same model as the one developed in Section [XTj except that an error term. 
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vl, is added to the output function. The errors {f^}^!^^ are assumed to be a sequence 
of independent normally distributed random variables with variance S. In this model 
structure they encompass measurement errors, as well as model specification errors and 
system stochasticity. When considering the whole population the main differences lie in 
K^, whose coefficients are 

4,r^N{0,ul,), (27) 
V\2-^N{0,iol,), (28) 

(29) 
(30) 

vUr^NiO^ul), (31) 
vl,r^N{0,u;l), (32) 
vUr^NiO,ul,), (33) 

(34) 



7^^^~iV(0,u;JJ. (35) 

The group index variable Gq (and G^^) equals one if individual i belongs to the control 
group (diabetic group) and zero otherwise. The fixed effects parameters k^^, kl2, k{^, kl-^^, 
^43, and Qi are the same for all individuals in group j and represent the parameters for a 
typical individual in respective group (where j is C or D). The parameters are assumed 
to have a lognormal distribution in the population to avoid negative parameter values. 

To sum up, the unknown population parameters that may be estimated in the ODE 
based NLME framework are: the structural population fixed effects parameters k^^, /c-Jg, 
k{^, kli, kl^, and Q{ (for j = C,D), the variances cugi, 1^125 1^135 1^31, 1^43, and Uq^ of 
the interindividual variation, and the variance S of the measurement error term. Observe 
that we assume that the interindividual and the intraindividual variations are the same 
for the two groups in this model. 

The results from the two stage approach is used as an initial guess of the structural 
parameters. The initial guess for 5* was obtained by first estimating 5* alone, with all the 
other parameters kept fixed to their initial values. 

3.2.1 Estimated parameters 

Table [2] shows the estimated parameters and standard errors (SEs) for the ODE model 
applied to the experimental data (the estimated parameters from the SDE model 
explained in the next section are also included in the table). Observe that the uncertainties 
of the cj-parameters are large compared to the estimated values. 





— (G^A;^! + G'^D^Q^ 


)exp(r7^i). 




= (G^/c^ + G^^k^2 


)exp(r7*2), 


'^L2 


= 0.01, 






— P 






= (G^fc^ + G^fcf^ 


)exp(r/l3). 


P 


= (G^fc^ + G^fc;^ 


)exp(r/^i). 


P 
^43 


= (G^A;^ + G^^k^^ 


)exp(7743). 


P 


= o.iki,, 





and in Q\, which is given by 

Ql = (G'cQf + G'^Qf)exp(r^^J, 



18 



In addition to estimating the parameters using all experimental data, we also estimated 
the parameters using a smaller set of data, in which the data points at 4, 8, 10, 15, 
and 45 minutes and 3, 6, and 8 hours were removed. Thus, at most eight samples per 
individual were used for the smaller set. The reason for estimating the parameters for a 
smaller data set was to test if population parameters can be reliably estimated even if 
the number of samples in the time series is low for each separate individual. It should 
be mentioned that the two stage approach was also applied to the smaller data set. In 
that case the spread of the estimated parameters was much larger between individuals 
compared to when all data were used and at least one parameter reached the predefined 
lower- or upper bounds for most individuals. 

In Table |2] and Figure H] it can be seen that using at most eight samples per individual 
gave no disjoint 85% confidence intervals compared to when all data were used (we used 
the standard errors of the estimated parameters to obtain confidence intervals for the 
NLME model). Most parameters are also very similar in general, but the estimates of 
^i2y ^3iy ^^cl wig deviates rather much for the two data sets and the overlapping confidence 
intervals are due to large uncertainties of the estimates. The estimated standard errors of 
the average parameters and the intraindividual noise parameter, S, were in general higher 
when only half the data set was used, indicating more uncertain parameter estimates. 

We were also interested in comparing the estimated population parameters from the two 





ODE - all data 


ODE - less data 


SDE - all data 




Controls 


Diabetics 


Controls 


Diabetics 


Controls 


Diabetics 


koi 


2.58 (0.13) 


2.08 (0.12) 


2.60 (0.13) 


2.05 (0.13) 


2.64 (0.14) 


2.06 (0.12) 


ki2 


1.75 (0.16) 


1.69 (0.20) 


1.51 (0.27) 


1.78 (0.22) 


1.61 (0.16) 


1.56 (0.20) 


ki3 


3.74 (0.54) 


5.47 (0.95) 


3.95 (0.57) 


4.93 (1.11) 


4.00 (0.50) 


4.82 (0.76) 




3.96 (0.39) 


3.40 (0.44) 


4.24 (0.44) 


2.88 (0.51) 


4.34 (0.48) 


3.32 (0.47) 


k43 


1.22 (0.17) 


1.38 (0.26) 


1.06 (0.23) 


1.54 (0.34) 


1.16 (0.14) 


1.27 (0.22) 


Qi 


297 (28) 


558 (59) 


296 (29) 


571 (64) 


290 (27) 


555 (59) 


^01 


0.044 (0.012) 


0.040 (0.011) 


0.044 (0.013) 


^12 


0.083 (0.034) 


0.015 (0.023) 


0.041 (0.034) 


^13 


0.258 (0.089) 


0.199 (0.098) 


0.073 (0.059) 


^31 


0.044 (0.040) 


0.024 (0.031) 


0.079 (0.043) 


^43 


0.296 (0.095) 


0.128 (0.071) 


0.203 (0.068) 




0.161 (0.040) 


0.173 (0.043) 


0.158 (0.041) 


s 


4.34- 10"^ 


0.38-10-^) 


4.89-10-^ 


0.74-10-^) 


6.98- 10-* (2.61 • 10-*) 


a 






3.09 (0.95) 


a 






0.222 (0.027) 



Table 2: Estimated parameters (and standard errors of the estimates) for the NLME model explained 
in equations ((27| - ((35|) . In the left panel the parameters of the ODE model (|25|l - (|26|) are presented 
when all experimental data were used. In the middle the estimated parameters are presented when 
the samples at 4, 8, 10, 15, and 45 minutes and 3, 6, and 8 hours were removed. The right panel 
consists of the corresponding estimated parameters, together with the parameters a and cr, when the 
SDE model dSS])-® was applied to all data. 
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stage approach with the results from the NLME approach when all data were used. As 
can be seen in Figure S] some estimated values differ between the two approaches, and 
the estimated medians fc^, and fcf^ have disjoint 85% confidence intervals. Moreover, 
all the interindividual variances are estimated to be smaller for the NLME model and 
Wqi, UJI2, ^13, and are smaller with statistical significance. 

Observe that the estimated population parameters for the two stage approach in Sec- 
tion [3XT] were trimmed (by removing estimates of parameters significantly different from 
most other individuals) before computing population averages and variations of the pa- 
rameters. If all estimated parameters are accounted for in that case the interindividual 
differences are even larger. Then the NLME approach predicts a significantly smaller 
interindividual variation also for /C43. 



Estimated medians 



Two stage, all data 
ODE model, all data 
ODE model, half data 
SDE model, all data 



uC hD 
^01 s-oi 



l,C l,D 



'"31 '^31 '"43 '"43 



Ql Ql 
100 100 



10| — 
8- 
6- 
4- 

2- 

— 



Estimated interindividual variability 
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Figure 4: 85% confidence intervals of the parameters that describe the medians of the estimated 
parameters and the interindividual variations for the three groups, respectively. Medians denoted by 
'x' represent the two stage approach (upper values for the intervals of lOcjj^j ^'^^ lOcjfg are 13.0 and 
12.7, respectively). The 'o':s represent the ODE based NLME approach applied to all data, the '-'is 
represent the ODE based NLME approach applied to the smaller set of data, and the 'n'is represent 
the SDE based NLME model applied to all data. 
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With the two stage approach we discovered that Qi differed between the groups. For the 
ODE based NLME model we see that there are significant differences for k^i and Qi. 

It is worth mentioning that the confidence intervals for the population parameters of the 
lognormal distributions in the NLME model may be calculated as in Section 13.1. by 
taking every individual into account. Doing so, the 85% intervals were in general more 
narrow compared to using the "standard error approach" and, in addition to km and 
Qi, also fci3 and have disjoint confidence intervals between the groups. Thus, when 
confidence intervals of the medians of the population parameters are computed by taking 
every individual into account, we conclude that both the leucine levels in the plasma and 
leucine kinetics differ between control individuals and individuals suffering from diabetes 
type 2. 



3.3 The SDE based NLME model 

We extend the ODE model by assuming that the tracee in compartment 1 is fluctuating 
around a steady state value. The state variable Q\ ^ is modelled with a mean-reverting 
Ornstein-Uhlenbeck process as discussed in Section 12.31 for a one-state compartmental 
model. The reason for testing the model under this assumption is that the amount of 
amino acid material is probably not exactly constant in the body, but varies over time. 



Daily rhythms of leucine levels have also earlier been reported (see, e.g.. lLavie and Lavie 



20061). 



The SDE model takes the form 

dQi, = a{Ql, - Q\,)dt + aQi.dWt, Q\iO) ~ N (^Ql„ ^^||^) , (36) 

dql = {K^qi)dt, q\0) = [ql„ 0, 0, 0]^, (37) 

Yi = log{^)+vl, vlr^N{0,S), (38) 

where Q\ q corresponds to Q\ in Section 13.21 Thus, compared to the ODE model, two 
new parameters, a and a, must be estimated. 



3.3.1 Estimated parameters 

When using the smaller set of data, the estimate of a converged to zero, and the SDE 
model was reduced to the ODE model. Therefore we only present the parameter estimates 
when the model was used on the complete data set. 

Overall, the SDE model produced similar parameter estimates as the ODE model (see 
Table [2] and Figure Sj). However, w^g) '^is^ ^ii: ^ deviate quite much, but only S have 
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overlapping confidence intervals. Interestingly, the parameters that describe interindivid- 
ual variability are, in general, smaller than for the ODE model. 

Like the ODE based NLME model, the SDE model predicts group differences for /cqi and 
Qi when the "standard error" approach was used to compute confidence intervals and for 
hoi, ki3, k^i, and Qi when all individual estimates were used. 

We observed strong correlations between the parameters describing the intraindividual 
variations of the model. The estimated correlations were 0.79 between a and a, —0.40 
between a and S, and —0.69 between a and S. 



3.4 Comparison between ODE and SDE results 



Both the ODE- and SDE based NLME models gave good fits to data (see the upper 
plots in Figure Ej) and small correlations between residuals (when SDEs are used the 
term residual refers to the observed data point subtracted by the one step prediction 
obtained from the extended Kalman filter at the given time point). The difference between 
the two approaches is that the SDE model predicts a varying leucine level. Let zl = 
Qi t/Qi be the normalized process described by dzl = a(l — zl)dt + adWt- Then zl has 
variance o"^/(2a), thus the leucine level is predicted to vary with a standard deviation of 
100 X a/ \/2a percent of the predicted average value. With the estimated values of a and 
a we get a standard deviation of 0.089. 

In Figure E] we visualize how the level of free leucine in the body changes in the extended 
Kalman filter for a single individual. First the model predicts the level to be at the mean 
level. Then an update is performed at the first measurement time point and the model 
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Figure 5: The residuals for all 34 individuals at the sixteen sampled time points. The lines correspond 
to mean values. In the lower plots the reduced models described in Section [3. 4. H are used. 
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makes a new prediction at the next time point, given the previous update. The procedure 
is then iterated until the last time point is reached. Observe that we at the x-axis have the 
index of each time point, not the actual times. In the plot the leucine level is normalized 
as described above. 



In the output functions, and (1551) . it is assumed that the errors vl are Gaussian and 
independent. The first assumption has been tested by applying the Lilliefors test function 
lillietest in MATLAB to the residuals. With significance level 0.05 an assumption about 
normality can be rejected for four individuals for the ODE model and one individual 
for the SDE model (in these individuals we observed one or two residuals significantly 
different from the other, indicating outlying observations). 

One way to investigate the independence assumption is by means of the aut ocorrelations 
between the residuals (one step pr ediction errors) . That was also done in lM0ller et al. 



( 120101 ) and lOvergaard et al.l ( 120071 ) . where they concluded that the residual correlations 
decreased for their proposed SDE models. 

However, the autocorrelation of the residuals does not give any information about corre- 
lations between specific sample points, but is more a general representation of the correla- 
tions for specific lags. When analyzing multiple individuals sampled at the same times we 
can instead compute the sample correlations of the residuals for specific time points (the 

correlation between residuals at time tj and tk is calculated as rjk = Xlili ^''\N-i)s sk''^ ^ 
where is the mean of all the residuals at time t^, and Sk is the sample standard devia- 
tion). This gives a more detailed representation of the residual correlation structure. For 



Normalized leucine level in the EKF for individual 1 in the control group 
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Figure 6: The normalized free leucine level in the EKF for an individual in the control group. The 
dashed lines indicate 1 ± 0.089, i.e., the mean plus/minus the standard deviation of the process at 
each time point. For this individual there is a missing data point at 8 hours, i.e., the data are taken at 
2, 4, 6, 8, 10, 12, 15, 20, 30, and 45 minutes and 1, 2, 3, 4, and 6 hours. Observe that when the 
interval between measurements increase, the jumps away from the mean becomes longer, indicating an 
increase in the uncertainty of the state. 
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the NLME models presented in Sections 13.21 and [ 3.31 the correlations between the residuals 
are small. A correlation matrix for the residuals from the ODE model is presented in the 
left hand plot in Figure [3 



3.4.1 Reduced ODE- and SDE based models 

If an ODE based model is less accurate it may produce subsequent over- or underpre- 
dictions of the data, thus producing correlated residuals. A method to reduce these 
correlations may be to introduce system noise by means of SDEs in the model. Since 
the ODE model described above seems to produce good fits to data we have investigated 
this for a reduced version of the leucine model where and c^is were fixed to zero (the 
reason to choose these parameters is that a sensitivity analysis indicate that ki2 and fcis 
have the least influence on the objective function). We call the new models the reduced 
ODE- and SDE models, respectively. 

The sample correlation matrices for the residuals of the reduced models are also presented 
in Figure [3 We observe that the residuals are less correlated for the SDE model. The 
most striking difference between the two models is that the correlations between adjacent 
residuals are in general positive for the ODE model but not for the SDE model. The 
conclusion to be drawn from this is that the SDE model well captures the involved mech- 
anisms and that the remaining noise and uncertainty can be described well by normally 
distributed random variables and Wiener processes. 



Unreduced SDE Model Reduced ODE Model Reduced SDE Model 




-0.5 0.5 1 -0.5 0.5 1 -0.5 O 0.5 1 



Figure 7: Correlation matrices for the residuals for the ODE based NLME model where all parameters 
vary between individuals, and the reduced ODE- and SDE models respectively. Only the correlations 
significantly different from zero [p < 0.05) are highlighted and correlations with p > 0.05 are set to 
zero. The p-values are calculated by testing the hypothesis that the correlation r^^ is zero against the 
alternative rj^. 7^ as explained in lLarsen and Marxl ( 200ll ). 
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4 Discussion 



In this article population models for leucine kinetics in blood plasma have been proposed. 
The models are based on differential equations to account for the dynamic behaviour of 
leucine molecules. The first models are based on ordinary differential equations where the 
tracee material is assumed to be in a steady state. The NLME model is then extended 
so that the amount of tracee material is modelled by an Ornstein-Uhlenbeck process to 
account for fluctuating tracee material. 

The estimated population average parameters did in general not change much when SDEs 
were used instead of ODEs in the NLME settings, but the interindividual variabilities were 
in general estimated to be smaller for the SDE model. In systems where the tracee level 
fluctuates much around a steady state during the time of the experiments, we believe 
that errors in the estimated physiological parameters may occur if the fluctuations are 
not accounted for in the model. We have performed simulation studies for the one- 
compartment SDE model (jl])-® that supports this hypothesis. Especially the parameters 
that describe the interindividual variability tend to be overestimated when an ODE based 
model is used for data generated from the SDE model. 

For the constrained model {uu = wis = 0) the residuals were correlated when ODEs were 
used. Correlations between the residuals indicate that the assumption that the errors 
Vk and Vi are independent when 7^ / is incorrect. For the constrained SDE model the 
correlations were smaller. This explains how SDE may be used as a tool to handle non- 
explicitly modelled effects in the system. For the present leucine kinetics model there 
is no obvious reason not to use the unconstrained model, but it should be mentioned 
that this model requires two more parameters to be estimated for each individual and 
the computational time needed to estimate the parameters with the unconstrained model 
increased significantly. For larger models it may be computationally too demanding to 
estimate all the parameters for each individual. 

Other S DE models rn a y be u sed to describe the system. One approach, that was recently 



used by lM0ller et al.l (120101 ). is to add a zero- mean Ornstein-Uhlenbeck process to the 
model output to allow for a more flexible model output structure. We have tested that 
approach and the results indicate that such a model may equally well be used to describe 
the system. However, that model does not have a direct physiological description and we 
decided to only present the results for the variable tracee model. 

The parameters of the ODE based NLME model were also estimated using half the data 
for each individual, and the results were similar to the ones obtained using all data. This 
indicates that the number of samples needed per individual does not necessarily have to 
be large for the ODE based NLME model. This may be of importance in experimental 
planning since taking many repeated measurements may be expensive or troublesome. 
The SDE model could not be applied when the smaller data set was used, since a converged 
to zero in that case, resulting in the ODE model. Thus, when the EKF is applied to the 
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SDE based leucin e kinetics model, t he nu mber of data points per individuals needs to be 
sufficiently large. lOvergaard et al.l (120051 ) show that, when the number of individuals is 
large, only three samples per individual are enough to estimate parameters for a simple 
one-compartment SDE model. Further studies are needed to investigate the number of 
individuals and data points per individual that is needed to estimate the parameters for 
our model. Of particular interest is to investigate the model performance for only a few 
data points per individual. To capture the dynamical behaviour, it is then important that 
the measurements are not taken at the same time points for each individual. 

Differences were observed when the estimated parameters from the two stage approach 
were compared to the results from the NLME models. Population median estimates dif- 
fered for some parameters but especially the interindividual variations were estimated 
to be much smaller for the NLME models. It has earlier been observed that, for some 
mo dels, the two stage app roach produce biased average parameters (see, e.g., the work 
by lSheiner and Beall 119831 ). That the interindividual variation is overestimated with the 
tw o stage appro a ch ha s earlier been obse r ved in a la rg e number of publica tions, e.g., 
in lOlofsen et aP tQ04 \: ISheiner and Beall f ll980l . Il983h : ISteimer et aJ f ll984h . One ex- 



planation why the two stage approach seems to produce biased estimates of population 
parameters is that the variance of the av erage parameters dep end both on intra- and in- 
terindividual variations. This is shown in lOlofsen et al.l (12004 ) for a simple linear model. 
NLME models separate these two kinds of variations. 

We have investigated the model parameters for two groups of individuals, one group 
of diabetic individuals and one control group. The estimated average behaviour of the 
parameters differs between the groups, and especially Qi is much larger for the diabetic 
group. Larger levels of free leucine for diabetic individuals have been reported earlier and 
the amount of leucine is associated with the level of insu lin treatment and magnitude of 
hyperglycemia for the individuals ( iGougeon et al.l . l2008l ) . Moreover leucine degradation 
and protein uptake of free leucine (/cqi and k^i) is smaller for the diabetic group and leucine 
production due to increased pro teolysis is large r . The se findings are also in agreement to 
what has been reported earlier ( Gougeon et al.l . 2008). 



There are also other physiological differences between the groups (age, body weight, level 
of obesity, etc.) that may influence the parameters and a thorough investigation how such 
covariates associates with model parameters is a topic of further research. Such knowledge 
is important when NLME results are used in predictive studies and when parameters 
are estimated for single individuals, especially when data are very sparse. The NLME 
framework can then be used by simply maximizing the individual loglikelihood (|T2|) with 
all the population parameters fixed to results from earlier population studies. Note that 
predictions of individual parameters can then be performed for very few data points (even 
a single data point) but the predictive performance obviously increase for individuals that 
have been densely sampled. 



To analyze the identifiability of the model parameters we have performed studies where 
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data were simulated for a given set of parameters. Given the simulated data we estimated 
the parameters and they were in general similar to the true parameters, except that a 
and cr were somewhat overestimated. 

It would be interesting to apply other filters, which arc adapted to nonlinear models. 
However, in this particular model it is only the output function that is nonlinear, and we 
believe that the EKF is a good choice. 

The use of nonlinear mixed effects models were originally used in the area of pharmacoki- 
netics and pharmacodynamics. Recently the number of applications has been growing 
and the present paper is a first step towards larger NLME models in the metabolic arena, 
focusing on studying lipoprotein kinetics. Using SDEs may be a good way to account for 
model approximations and/or true variations of model state variables and parameters. 
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